#include "Matrix.h"
#include "Vector.h"
#include "DeterminantCalculator.h"
#include "LinearSolver.h"

Matrix::Matrix ( size_t rowsCount, size_t columnsCount)
{
	rows = new Vector*[rowsCount];

	for (size_t i = 0; i < rowsCount; i++)
	{
		rows[i] = new Vector (columnsCount);
	}

	this->rowsCount = rowsCount;
	this->columnsCount = columnsCount;
}

Matrix::Matrix ( size_t rowsCount, size_t columnsCount, const Vector *vector)
{
	if (NULL == vector)
	{
		throw std::invalid_argument("Pointer to vector array cannot be NULL");
	}
	
	rows = new Vector*[rowsCount];

	for (size_t i = 0; i < rowsCount; i++)
	{
		rows[i] = new Vector (columnsCount, vector[i]);
	}
	
	this->rowsCount = rowsCount;
	this->columnsCount = columnsCount;
}

Matrix::Matrix ( size_t rowsCount, size_t columnsCount, const double *data )
{
	if(NULL == data)
	{
		throw std::invalid_argument("Pointer to data array cannot be NULL");
	}
	rows = new Vector*[rowsCount];

	for(size_t i = 0; i < rowsCount; i++)
	{
		rows[i] = new Vector(columnsCount, data + columnsCount * i);
	}

	this->rowsCount = rowsCount;
	this->columnsCount = columnsCount;
}

Matrix::Matrix (const Matrix &another)
{
	rows = new Vector* [another.rowsCount];
	
	for (size_t i = 0; i < another.rowsCount; i++)
	{
		rows[i] = new Vector(*(another.rows[i]));
	}
	
	rowsCount = another.rowsCount;
	columnsCount = another.columnsCount;
}

Matrix::Matrix ( size_t rowsCount, size_t columnsCount, const Matrix& another)
{
	rows = new Vector* [rowsCount];

	for (size_t i = 0; i < another.rowsCount; i++)
	{
		rows[i] = new Vector(columnsCount, *another.rows[i]);
	}
	
	for (size_t i = another.rowsCount; i < rowsCount; i++)
	{
		rows[i] = new Vector(columnsCount);
	}

	this->rowsCount = rowsCount;
	this->columnsCount = columnsCount;
}

size_t Matrix::getColumnsCount() const
{
	return columnsCount;
}

size_t Matrix::getRowsCount() const
{
	return rowsCount;
}

Vector Matrix::column (size_t columnIndex) const
{
	if ( columnIndex > columnsCount - 1)
	{
		throw std::out_of_range("Column index cannot be more then matrix width");
	}

	Vector vector(rowsCount);

	for (size_t i = 0; i < rowsCount; i++)
	{
		vector[i] = (*rows[i])[columnIndex];
	}

	return vector;
}

Matrix::~Matrix ()
{
	for (size_t i = 0; i < rowsCount; i++)
	{
		delete rows[i];
	}

	delete[] rows;
}

const Vector& Matrix::operator [] ( size_t index) const
{
	if ( index > (rowsCount - 1) && index < 0 )
	{
		throw std::out_of_range("");
	}
	return *(rows[index]);
}

Vector& Matrix::operator [] ( size_t index)
{
	if ( index > (rowsCount - 1) && index < 0 )
	{
		throw std::out_of_range("");
	}
	return *(rows[index]);
}

Matrix& Matrix::operator = (const Matrix& another)
{
	if (this == &another)
	{
		return *this;
	}

	for (size_t i = 0; i < rowsCount; i++)
	{
		delete rows[i];
	}

	delete[] rows;

	rows = new Vector* [another.rowsCount];
	
	for (size_t i = 0; i < another.rowsCount; i++)
	{
		rows[i] = new Vector(*(another.rows[i]));
	}
	
	rowsCount = another.rowsCount;
	columnsCount = another.columnsCount;

	return *this;
}

std::ostream& operator << (std::ostream& outputStream, const Matrix& matrix)
{
	outputStream << "{";
	
	for( size_t i = 0; i < matrix.rowsCount - 1; i++ ) {
		outputStream << *(matrix.rows[i]) << ", ";
	}
	if (matrix.rowsCount > 0)
	{
		outputStream << *(matrix.rows[matrix.rowsCount - 1]);
	}
	outputStream << "}";
	return outputStream;
}

std::istream& operator >> (std::istream& inputStream, Matrix& matrix)
{
	if( 0 == matrix.rowsCount ) {
		return inputStream;
	}

	Matrix matrixToRead (matrix.rowsCount, matrix.columnsCount);

	char  matrixSeparator;
	inputStream >> matrixSeparator;
	if( inputStream.fail() || matrixSeparator != '{' ) {
		throw std::invalid_argument ("Matrix begin is not {");
	}

	char delimiter;
	for(size_t i = 0; i < matrixToRead.rowsCount - 1; i++ ) {
		inputStream >> *(matrixToRead.rows[i]);
		inputStream >> delimiter;
		if( inputStream.fail() || delimiter != ',' ) {
			throw std::invalid_argument ("Matrix components must have a delimiter ,");
		}
	}
	inputStream >> *(matrixToRead.rows[matrixToRead.rowsCount - 1]);
	inputStream >> matrixSeparator;
	if( inputStream.fail() || matrixSeparator != '}' ) {
		throw std::invalid_argument ("Matrix end is not }");
	}
	
	matrix = matrixToRead;
	return inputStream;
}

Vector operator* (const Matrix& matrix, const Vector& vector)
{
	if ( matrix.columnsCount != vector.size() )
	{
		throw std::invalid_argument("Matrix width must be equal to vector size");
	}
	Vector result ( matrix.rowsCount );

	for ( size_t i = 0; i < matrix.rowsCount; i++)
	{
		result[i] = *(matrix.rows[i]) * vector;
	}
	return result;
}

bool Matrix::isValidRow(size_t row) const {
	return ( row < rowsCount && row >= 0 );
}

bool Matrix::isValidColumn(size_t column) const {
	return ( column < columnsCount && column >= 0 );
}

void Matrix::exchangeRows(size_t firstRow, size_t secondRow) {
	if(!isValidRow(firstRow) || !isValidRow(secondRow) ) {
		throw new std::invalid_argument("Row index is incorrect");
	}

	Vector *tmp = rows[firstRow];
	rows[firstRow] = rows[secondRow];
	rows[secondRow] = tmp;
}

void Matrix::exchangeColumns(size_t firstColumn, size_t secondColumn) {
	if(!isValidColumn(firstColumn) || !isValidColumn(secondColumn) ) {
		throw new std::invalid_argument("Column index is incorrect");
	}
	for (size_t i = 0; i < rowsCount; i++) {
		rows[i]->exchange(firstColumn, secondColumn);
	}
}

Matrix Matrix::transpose() const
{
	Matrix result(rowsCount, columnsCount);

	for(size_t  i = 0; i < rowsCount; i++)
	{
		for (size_t j = 0; j < columnsCount; j++)
		{
			result[i][j] = (*rows[j])[i];
		}
	}

	return result;
}

double Matrix::determinant() const
{
	if( rowsCount != columnsCount ) {
		throw std::invalid_argument("Determinate calculation allowed only for square matrix");
	}
	if( 0 == rowsCount ) {
		return 0.0;
	}
	return DeterminantCalculator::calculate(*this);
}

Vector Matrix::solveLinearEquation (const Vector& freeTermsVector)
{
	if( rowsCount != columnsCount ) {
		throw std::invalid_argument("We can solve linear equation for quad matrix only");
	}
	if( freeTermsVector.size() != columnsCount ) {
		throw std::invalid_argument("Vector size must be equal to matrix columns count");
	}
	return LinearSolver::solve (*this, freeTermsVector);
}


Matrix operator* (double constant, const Matrix& another)
{
	Matrix result (another);

	for ( size_t i = 0; i < result.rowsCount; i++ )
	{
		for (size_t j = 0; j < result.columnsCount; j++)
		{
			result[i][j] *= constant;
		}

	}
	return result;
}

Matrix operator* (const Matrix& another, double constant)
{
	return operator*(constant, another);
}

Vector operator*(const Vector& vector, const Matrix& matrix)
	{
		
		Vector result(matrix.columnsCount);

		for(size_t i = 0; i < matrix.columnsCount; i++)
		{
			result[i] += vector * matrix.column(i);
		}

		return result;
	}

Matrix operator* (const Matrix& first, const Matrix& second)
{
	
	Matrix result(first.rowsCount, second.columnsCount);

		for(size_t i = 0; i < first.rowsCount; i++)
		{
			result[i] = first[i] * second;
		}

		return result;
}

Matrix Matrix::inverse () const
{
	size_t  size = rowsCount;
	Matrix tmpMatrix(size, 2 * size);
	Matrix I(size, size);
	for(size_t i = 0; i < size; i++)
	{
		for(size_t j = 0; j < size; j++)
		{
			tmpMatrix[i][j] = (*rows[i])[j];
		}
		tmpMatrix[i][i + size] = 1;
	}
	for(size_t i = 0; i < size - 1; i++)
	{
		int pivot = i;
		for(size_t j = i + 1; j < size; j++)
		{
			if(abs(tmpMatrix[j][i]) > abs(tmpMatrix[pivot][i]))
			{
				pivot = j;
			}
		}
		if(pivot != i)
		{
			Vector tmpVector = tmpMatrix[i];
			tmpMatrix[i] = tmpMatrix[pivot];
			tmpMatrix[pivot] = tmpVector;
		}
		for (size_t j = i + 1; j < size; j++)
		{
			if(0 == tmpMatrix[i][i])
			{
				throw std::invalid_argument("Matrix is singular");
			}
			double factor = tmpMatrix[j][i] / tmpMatrix[i][i];
			for(size_t k = 0; k < 2 * size; k++)
			{
				tmpMatrix[j][k] -= tmpMatrix[i][k] * factor;
			}
		}
	}
	if(0 == tmpMatrix[size - 1][size - 1])
	{
		throw std::invalid_argument("Matrix is singular");
	}
	for(size_t column = 0; column < size; column++)
	{
		for(size_t i = size - 1; i >= 0; i--)
		{
			double sum = 0;
			for(size_t j = i + 1; j < size; j++)
			{
				sum += tmpMatrix[i][j] * I[j][column];
			}
			I[i][column] = (tmpMatrix[i][column + size] - sum) / tmpMatrix[i][i];
		}
	}
	return I;
}